Filling the holes: Evolving excised binary black hole initial data with puncture 

techniques 
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We follow the inspiral and merger of equal-mass black holes (BHs) by the moving puncture tech- 
nique and demonstrate that both the exterior solution and the asymptotic gravitational waveforms 
are unchanged when the initial interior solution is replaced by constraint-violating "junk" initial 
data. We apply this result to evolve conformal thin-sandwich (CTS) binary BH initial data by 
filling their excised interiors with arbitrary, but smooth, initial data and evolving with standard 
puncture gauge choices. The waveforms generated for both puncture and fiUed-CTS initial data are 
remarkably similar, and there are only minor differences between irrotational and corotational CTS 
BIf binaries. Even the interior solutions appear to evolve to the same constraint-satisfying solution 
at late times, independent of the initial data. 



I. INTRODUCTION 

Numerical relativity has made dramatic progress over 
the past two years in solving the binary black hole (BBH) 
problem via the "moving puncture" technique P, H, [1| . 
To construct puncture initial data, the spatial metric is 
taken to be conformally flat and the conformal factor is 
split into an algebraic singular piece and a numerically 
determined regular contribution f3, 'Bj . These initial data 
are then evolved via the BSSN [6, 7] formulation of Ein- 
stein's equations, without special treatment of the initial 
(coordinate) singularities at the punctures. Coupled to 
the BSSN field evolution equations are a set of gauge 
evolution equations, similar to those discussed in [8|, |9| 
and successfully evolved by [H, that are critical for 
long-term stability. 

One possible reservation is that puncture initial data 
are thought to contain more inherent eccentricity (lo| 
and, possibly, more initial spurious gravitational radia- 
tion than those constructed using the conformal thin- 
sandwich (CTS) method. CTS initial data are typically 
constructed with equilibrium boundary conditions im- 
posed on an excised surface defining the BH horizon 
[111 [H, [13], and it has been suggested that the lack 
of data in the BH interior may pose a significant prob- 
lem for evolving CTS initial data with puncture tech- 
niques (see, e.g., Thus, excision-based initial data, 
for both orbiting binaries [la, [lH and head-on collisions 
[itI . [TH , have heretofore been evolved via dynamical ex- 
cision, rather than puncture techniques. We will demon- 
strate here that there are no difficulties evolving BH inte- 
riors artificially filled with smoothly extrapolated "junk" 
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(i.e., constraint- violating) initial data via puncture tech- 
niques: the puncture gauge conditions quickly drive the 
evolution toward the usual late-time puncture interior 
solutions, independent of the initial interior data. 

This demonstration follows from several previous stud- 
ies of BHs in puncture gauges. In [l^ [13, [Mj, it was 
shown that the conformal factor ip evolves from a 1/r 
dependence at small r to l/xA" at late times, indicating 
that the spatial slices terminate on a limiting surface of 
finite areal radius. It was argued in that this effect is 
really a consequence of "excision via under-resolution" , 
namely the fact that the puncture gauge conditions make 
numerical grid points move away from the puncture very 
rapidly. In a Kruskal diagram, where the puncture is 
represented by the asymptotic end of region III, all grid 
points quickly move from region III into region II, leav- 
ing few or no grid points behind in the "left" part of the 
diagram (see Fig. 2 in (2^). 

In a recent paper we demonstrated that one can 
replace the solution inside the horizon of a single station- 
ary puncture BH by "junk data" . We find that this junk 
does not affect the BH exterior and, more surprisingly, 
even in the black hole interior the BH settles down to the 
same constraint-satisfying solution, independently of the 
initial data. This can again be motivated with the help 
of a Kruskal diagram. We impose junk on the t = slice 
in the black hole interior, i.e. in the "left" part of the dia- 
gram. As demonstrated by [2^, the puncture gauge con- 
ditions make all grid points propagate toward the right at 
extremely rapid, "superluminal" speeds. Assuming that 
constraint violations travel with much smaller speeds, the 
grid points are likely to leave the junk's future domain of 
dependence eventually, and enter a region that is affected 
by constraint-satisfying initial data only. 

Here, we extend our work to the BBH case. We first 
demonstrate numerically that replacing interior puncture 
initial data by "junk" data does not alter the exterior 
evolution. Given this finding, we then show that the 
excised CTS initial data generated by Cook and Pfeiffer 
([HI; hereafter CP) can be evolved using puncture gauges 
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without loss of accuracy by filling in the interior with ar- 
bitrary, but smooth, initial data. Rather than producing 
highly accurate waveforms, the goal of this paper is to 
establish a point of principle about initial data for dy- 
namical evolutions. 



II. INITIAL DATA 

Throughout this paper we cast the spacetime metric 
gab into its 3+1 ADM form 

= _(q-2 - PiP')dt^ + 2l3idtdx' + i^'^^ydx^dx^ , (1) 

with lapse function a, shift vector /3% and spatial metric 
7y = where jij is the conformally related metric 

and the conformal factor ip is chosen so that det(7) — 1 
in Cartesian coordinates. Both puncture and CP spatial 
initial data are conformally flat at t — 0, i.e., jij = rjij. 

In generating BBH puncture initial data 0, Q, the 
momentum constraint reduces to an algebraic expression 
and the Hamiltonian constraint to a nonlinear elliptic 
equation governing the nonsingular part of tp, which we 
solve via the Lorene multidomain spectral methods li- 
braries [1^. These constraint equations depend on the 
initial linear momenta, coordinate separation, spins, and 
puncture masses. We obtain these initial parameters 
from the quasiequilibrium, circular sequence of ^25^, first 
treating the case of nonspinning, equal-mass BHs. 

To better understand the role played by the BH in- 
teriors in puncture evolutions, we replace the puncture 
interiors ("P") with "junk" initial data. Specifically, wc 
replace the singular term of ip with some other smooth 
function inside a radius r j and leave the nonsingular part 
untouched, as we did in [231 for stationary BHs. 

The singular part of V' is i/'s = where Mi 

is the puncture mass of BH i, and is the distance from 
the origin to BH i. We construct smooth junk ("Psj") by 
replacing the term Mi/2ri by an even, fourth-order poly- 
nomial inside rj /Mi — 0.4 (in these isotropic coordinates, 
the apparent horizon radius is located at tah ~ Q.hMi 
initially) , chosen such that the resulting conformal factor 
is continuous and twice differentiable everywhere. Flat 
junk ("Ppj") is defined by replacing Mijlvi with a con- 
stant value 2.0 within rj/Mi = 0.25. This second choice 
produces a conformal factor that is continuous, but not 
differentiable everywhere. Clearly, the the Hamiltonian 
constraint is not satisfied within the junk regions of our 
grid initially. 

The corotating CP initial data ("CPc") are available 
from 26], whereas the irrotational CP data ("CPi") were 
graciously supplied by Harald Pfeiffer. Unlike puncture 
initial data, which are defined everywhere in space, CP 
initial data are produced by imposing equilibrium bound- 
ary conditions on excision surfaces - chosen to be coor- 
dinate spheres ~ that define the BH's apparent horizons. 
This procedure produces valid data only in the BH ex- 
terior. To aide in some numerical simulations, the CP 
initial data were provided with data extrapolated from 



TABLE L Summary of parameters for our simulations. For 
initial data, "CP" refers to CTS Cook-Pfeiffer data, smoothly 
extrapolated inside the excised regions. M is the total ADM 
mass of the binary system, ro and Jo the initial separation 
and total angular momentum, and Mirr the irreducible mass 
of each BH. A_B and AJ are the total radiated energy and 
angular momentum during the course of the evolution in the 
dominant I = 2, m = ±2 modes, evaluated at extraction radii 
of 17.8M. 



Name Init. Data 


MQ. ra/M .h/M'^ M^rr/M 


AE/AI AJ/Jo 


P No Junk 
Psj Sm. Junk 
Ppj Flat Junk 


0.0824 2.185 0.807 0.476 
0.0824 2.185 0.807 0.476 
0.0824 2.185 0.807 0.476 


0.026 0.16 
0.027 0.17 
0.026 0.16 


CPi Irrot. CP 
CPc Corot. CP 


0.0824 2.237 0.823 0.447 
0.0817 2.222 0.877 0.444 


0.029 0.19 
0.028 0.19 



the BH exterior to a region outside O.TSr-AH- To fill the 
remaining excised region inside O.TSrAH? we extrapolate 
all field values radially from O.SrAH to the BH center, 
matching smoothly to fourth-order polynomials with a 
fixed, finite value at r = 0, so that the resulting field 
variables are everywhere continuous and twice differen- 
tiable, similar to the Psj puncture case. Evolving with 
the CTS gauge variables a and results in eccentric 
BH coordinate trajectories. We can remove this gauge- 
dependent eccentricity considerably by choosing instead 
the same initial lapse and shift as in the puncture cases, 
a = -0"^ and = 0, evolving with Eqs. © - dS]) below. 

In Table U we summarize the parameters for each 
of our numerical simulations. Both irrotational and 
corotating CP data are evolved from an initial separa- 
tion To ~ 2.2M, which is beyond the quasiequilibrium 
ISCO radius, and comparable to the initial configura- 
tions evolved in [J, 0] ■ Here and below M is the binary's 
total ADM mass. For comparison, we construct puncture 
initial data with the same orbital frequency as the irro- 
tational CP data {Mil = 0.08241), using the equal-mass 
irrotational fitting formulae of 25'| to determine our ini- 
tial parameters. All runs were performed on 420^^ x 210 
spatial grids, assuming equatorial symmetry. We employ 
multiple-transition fisheye coordinates [l| to expand the 
physical extent of our grid, as we did in (23l |. Our grid 
includes three 1:4 fisheye transitions at r/M = 3.89, 5.22, 
and 6.55, all with transition width s/M = 0.417. Resolu- 
tion in the innermost region is set to M/24, so the outer 
boundary falls at a physical distance of « 170M. 



III. DYNAMICAL CALCULATIONS 

Our evolutions are performed using the same Cactus- 
based [23 GRMHD code described in [H, but with 
the MHD matter and E&M field sectors disabled. Our 
originally second-order Iterative Crank- Nicholson BSSN- 
based code has been upgraded to fourth-order accurate 
spatial differencing with upwind differencing for advec- 



3 



tive terms, but we still use second-order differencing for 
the time evolution. 

We use standard puncture gauge evolution equations, 
specifically a 1 + log lapse and a "non-shifting shift" : 

dta = P^dja-2aK (2) 
dtl3' = (3/4)5^ dtB' = dtV - vB' (3) 

where we set 77 — 0.5 /M in all of our runs. 




so 40 60 



t/M 

FIG. 1: Binary separation versus time and patlis traced in 
coordinate space (inset) for runs P (thin solid), Pfj (thin 
dashed), Psj (thin dotted), CPi (thick solid), and CPc (thick 
dashed) . In the inset, we show the motion of one of the BHs in 
the equatorial plane, with points representing positions every 
lOM for runs P (open squares) and CPi (closed squares). 
While TT-symmetry is not enforced, it holds true numerically 
during our evolutions. 

In Fig. [U we show the coordinate trajectories of one 
of the punctures, which obey the relation dx^/dt — 
— /3*(a;')[]|. For runs with non-singular junk data in the 
interior, we locate the initial BH position at r = 0, coin- 
cident with the standard puncture case. While the punc- 
ture and CP runs trace out slightly different tracks, this 
comparison is gauge-dependent. To evaluate our runs in 
an asymptotically gauge-invariant manner, we compute 
the gravitational waveforms through the Weyl scalar ^^4 
via the Cactus PsiKadelia thorn 28], which we have 
modified to make the computation of spatial derivatives 
accurate to fourth order. We have tested our modified 
PsiKadelia thorn with a full suite of analytic linearized 
wave tests. The total energy and angular momentum 
lost to gravitational radiation are computed by integrat- 
ing Eqs. (5.1) and (5.2) in [2§|, keeping only the s = — 2 
spin- weighted, / = 2, m = ±2 contributions from 
to reduce numerical noise. We have performed simula- 
tions at various resolutions using the same puncture ini- 
tial data as in ^ , finding that our waveform agrees with 
Fig. 3 of the preprint version of [1] to high accuracy, and 
converges to second order, as we show in Fig. [3] for the 
CPi case. 

In [131 , we demonstrated that even if single stationary 
punctures are constructed with different junk initial in- 
terior data, the ADM mass versus time and conformal 
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FIG. 2: Real part of the s = —2 spin-weighted, I — 2, m — 2 
component of the Weyl scalar ^4, evaluated on a sphere of 
radius r = 17. 8M. In the top panel we show results for punc- 
ture runs P (solid), Pfj (dashed), and Ps.j (dotted), which 
differ in the form of the initial junk we impose within each 
BH horizon. In the bottom panel, we compare run P with 
our two runs started from CP initial data, both irrotational 
(dashed) and corotating (dotted). We see extremely close 
agreement between puncture and irrotational CP results, and 
only minor differences between these and the corotating CP 
results. 
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FIG. 3: Waveforms for runs using the initial data of run 
CPi (top panel) evaluated at numerical resolutions of M/24 
(solid), M/20 (dashed), and M/16 (dotted). Results are time 
and phase shifted so that the moment and phase of the max- 
imum wave amplitudes coincide. In the bottom panel, we 
show the differences between the runs, properly rescaled to 
highlight second-order convergence. 



factor profiles at late times agree to high accuracy. Here, 
we demonstrate that the gravitational waveforms gener- 
ated by puncture BBHs are also largely independent of 
interior junk. In Fig. [21 we show the real Z = 2, m = 2- 
mode component of ■04, evaluated on a sphere of physical 
(as opposed to fisheye) radius r = 17. 8M, comparable to 
the extraction radius used in [H, 0] . In the top panel, we 
show a comparison between simulations P, Ppj, and Psj- 
We find that the first two cases are almost identical, even 
though the data Ppj are not differentiable initially in the 
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BH interior. There is a slight phase difference for run 
Pgj, since the differencing stencil in the inner edge of the 
exterior overlaps junk data in the interior that do not ini- 
tially satisfy the constraints. In the bottom panel of the 
figure, we compare results from simulation P with those 
from CP initial data (the CP curves are time and phase- 
shifted equally so that the time and phase at maximum 
amplitude match between runs P and CPi). We note that 
after an early phase of transient gravitational radiation 
propagation, the waveforms from the merger and ring- 
down of runs P and CPj match very closely, varying in 
amplitude by no more than 4% with nearly exact phase 
agreement. These relative differences are well within the 
relative differences in the initial data's masses and angu- 
lar momenta, as listed in Table HI 

The similarity between puncture and CP merger wave- 
forms may not be particularly surprising, since the gauge 
evolution equations inevitably drive a BH toward the uni- 
versal form described in To demonstrate this effect, 
in Fig. U we compare profiles of rilP' for runs P and CPj 
at t = and 20M, before coalescence. In both cases, 
the gauge conditions we choose drive the solution toward 
ip oc r~^^^, regardless of whether the initial data are sin- 
gular or expressly finite. For an isolated BH, it was shown 
in [l^ that the puncture r = evolves to a limiting slice 
of finite Schwarzschild radius, rs = 'ip'^r « 1.3M. We ob- 
serve the same asymptotic value for this function in all 
of our binary evolutions prior to coalescence. 
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FIG. 4: Radial profiles of the quantity ri/)^ as a function of 
coordinate radius, measured out from the BH center for runs 
P and CPi. The evolution inevitably drives ?/) toward the 
power-law index tp oc r~^^^ derived in [l^. For the curves 
showing profiles at t = 20M, we eliminate the two innermost 
points, whose values are inaccurate due to differencing across 
the puncture. Crosses denote the initial apparent horizons 
for each run. In the inset, we plot the value of r%lP linearly 
extrapolated from the third and fourth innermost data points 
as a function of time for the same runs. The filled circles show 
the asymptotic value 1.3M determined analytically by [l9l |. 

In Fig. [5l we compare the magnitude of the normal- 
ized constraint violations [s^l for different initial data 
sets, summing over spherical shells about each puncture 
spanning 0.2 < r/M < 0.3. These shells allow us to study 




FIG. 5: Violation of the Hamiltonian {7i — 0) and mo- 
mentum {A4^ — 0) constraints. We show the normalized x- 
component of the momentum constraint violation < > = 
\/'^{M^)'^ /Nmc (top panel), where Nmc is the L2-norm of 
A4^ given by Eq. 60 of ^3l> and a similar quantity for the 
Hamiltonian constraint (bottom panel). These are summed 
over spherical shells about each puncture 0.2 < r/AI < 0.3, 
which extend into the BH interiors. The constraint-violating 
initial data are driven to constraint-satisfying solutions at late 
times (several M). 



constraint violations propagating outward in the immedi- 
ate vicinity of the BHs, the exterior region most sensitive 
to the use of interior junk. We find a brief but transient 
increase in the magnitude of the violations at early times 
in all the puncture simulations, whereas the CP simula- 
tions begin with larger violations but are driven toward 
constraint-satisfying configurations quite rapidly. 



IV. DISCUSSION 

For solutions of Einstein's equations, the BH exterior 
must, by definition, be independent of the BH interior. 
As long as constraint violations propagate "causally", 
this must also be true if the BH interior violates the con- 
straints (see the subsequent discussion in 31]). Previous 
analyses [l^ [13, HT], [l^] also suggest that, with the "mov- 
ing puncture" gauges, any initial BH configuration would 
be driven to the asymptotic puncture solution found in 
. We demonstrate numerically that this is indeed the 
case and show that one can modify the data within the 
BH interior and produce valid waveforms. The success 
of the puncture method does not depend on the form of 
the initial data, or even on the existence of valid initial 
data, in the BH interior. 

More specifically, we find that evolutions of puncture 
initial data are insensitive to interior "junk" so long 
as the finite differencing stencils of exterior points do 
not overlap regions where constraints are violated. One 
may freely modify the interior fields to eliminate initially 
asymptotic singular behavior, or even invent initial data 
in the interior if they do not exist, so long as the solution 



5 



matches smoothly to the exterior so that all fields are 
twice difFerentiable and satisfy the constraints in regions 
covered by the differencing stencils of exterior points. 

As a consequence, excised initial data, in particular the 
CTS sets constructed in [ll], [l^, [3 , may be evolved us- 
ing puncture techniques. These data produce a different 
initial burst of spurious gravitational radiation, reflecting 
small differences in the "imperfection" of the initial data, 
but evolve similarly to punctures thereafter. This is true 
for both irrotational and corotating CTS data. The dif- 
ferences are small since the latter involve only marginally 
spinning BHs {S/M?.^ = 0.15). In fact, all our waveforms 
are very similar, and any relative differences in the wave 
signals are well within the relative differences of the phys- 
ical parameters describing the initial data. 

These results are also highly relevant for black hole- 
neutron star (BHNS) binary evolutions. To date, most 
published relativistic quasicquilibrium BHNS sequences 
are calculated using the CTS method and excision [32,. 
l33l . 1341 (the exception are the puncture data constructed 
in [3^[3g|). As demonstrated in this paper, these data 



can nevertheless be evolved with the moving puncture 
technique by appropriately filling the excised region. 

Finally, we note that the interior metric seems to evolve 
to the same asymptotic form found by ^19*1 , independently 
of the initial interior data. Consequently, the late-time 
interior solutions on the numerical grid satisfy the con- 
straints with increasing accuracy, even if the initial data 
consist of constraint- violating "junk" (see Fig. This 
feature presumably helps stabilize the integrations and 
explains the robustness of the puncture approach. 
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doin College. All simulations were performed on the 
NCSA abe cluster. 



[1] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo- 
chower, Phys. Rev. Lett. 96, 111101 (2006). 

[2] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006). 

[3] For an alternate method, see F. Pretorius, Phys. Rev. 
Lett. 95, 121101 (2005). 

[4] S. Brandt and B. Briigmann, Phys. Rev. Lett. 78, 3606 
(1997). 

[5] T. W. Baumgarte, Phys. Rev. D 62, 024018 (2000). 
[6] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 
(1995). 

[7] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 
024007 (1999). 

[8] C. Bona, J. Masso, E. Seidel, and J. Stela, Physical Re- 
view Letters 75, 600 (1995), arXiv:gr-qc/9412071. 
[9] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, 

D. PoUney, E. Seidel, and R. Takahashi, Phys. Rev. D 
67, 084023 (2003), arXiv:gr-qc/0206072. 

[10] E. Berti, S. Iyer, and C. M. Will, Phys. Rev. D 74, 
061503 (2006). 

[11] G. B. Cook and H. P. Pfeiffer, Phys. Rev. D 70, 104016 
(2004). 

[12] M. CaudiU, G. B. Cook, J. D. Grigsby, and H. P. Pfeiffer, 

Phys. Rev. D 74, 064011 (2006). 
[13] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom, 

G. Lovelace, and M. A. Scheel, Class. Quant. Grav. 24, 

59 (2007). 

[14] S. Husa, J. A. Gonzalez, M. Hannam, B. Bruegmann, 
and U. Sperhake, ArXiv e-prints (2007), 0706.0740. 

[15] J. G. Baker, M. Campanelli, F. Pretorius, and Y. Zlo- 
chower, Class. Quant. Grav. 24, 25 (2007). 

[16] J. Thornburg, P. Diener, D. PoUney, L. RezzoUa, 

E. Schnetter, E. Seidel, and R. Takahashi, ArXiv e-prints 
(2007), gr-qc/0701038. 

[17] U. Sperhake, ArXiv e-prints (2006), gr-qc/0606079. 
[18] U. Sperhake, B. Bruegmann, J. Gonzalez, M. Hannam, 



and S. Husa, ArXiv e-prints (2007), 0705.2035. 
[19] M. Hannam, S. Husa, D. PoUney, B. Briigmann, and 

N. O'Murchadha, ArXiv e-prints (2006), gr-qc/0606099. 
[20] M. Hannam, S. Husa, B. Briigmann, J. A. Gonzalez, 

U. Sperhake, and N. O. Murchadha, J. Phys. Conf. Series 

66, 2047 (2007). 
[21] T. W. Baumgarte and S. G. Naculich, Phys. Rev. D 75, 

067502 (2007). 
[22] J. D. Brown, ArXiv e-prints (2007), 0705.1359. 
[23] J. A. Faber, T. W. Baumgarte, Z. B. Etienne, S. L. 

Shapiro, and K. Taniguchi, Phys. Rev. D (2007), sub- 
mitted, 0708.2436. 
[24] http: //www . lorene . obspm.fr/. 

[25] W. Tichy and B. Briigmann, Phys. Rev. D 69, 024006 
(2004). 

[26] http: //www. black-holes . org/researchers3 .html. 
[27] http://www.cactuscode.org/. 

[28] We define the normalization of the tetrad to agree with 
that used in 15]. 

[29] J. Baker, M. Campanelli, C. O. Lousto, and R. Taka- 
hashi, Phys. Rev. D 65, 124012 (2002). 

[30] M. D. Duez, P. Marronetti, S. L. Shapiro, and T. W. 
Baumgarte, Phys. Rev. D 67, 024004 (2003). 

[31] D. Brown, O. Sarbach, E. Schnetter, M. Tiglio, P. Di- 
ener, I. Hawke, and D. PoUney, ArXiv e-prints (2007), 
0707.3101. 

[32] K. Taniguchi, T. W. Baumgarte, J. A. Faber, and S. L. 

Shapiro, Phys. Rev. D 74, 041502(R) (2006). 
[33] K. Taniguchi, T. W. Baumgarte, J. A. Faber, and S. L. 

Shapiro, Phys. Rev. D 75, 084005 (2007). 
[34] P. Grandclement, Phys. Rev. D 74, 124002 (2006). 
[35] M. Shibata and K. Uryu, Class. Quant. Grav. 24, 125 

(2007). 

[36] M. Shibata and K. Uryu, Phys. Rev. D 74, 121503 
(2006). 



